package fr.inserm.um915.bdibd.model;

import java.util.Set;
import java.util.SortedMap;
import java.util.TreeMap;
import java.util.TreeSet;

import org.apache.commons.math.MathException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.distribution.BinomialDistributionImpl;
import org.apache.commons.math.distribution.ContinuousDistribution;
import org.apache.commons.math.distribution.Distribution;
import org.apache.commons.math.distribution.IntegerDistribution;
import org.apache.commons.math.distribution.NormalDistributionImpl;

import com.sleepycat.je.Transaction;



/*
  
http://code.google.com/p/renjin/source/search?q=Distributions&origq=Distributions&btnG=Search+Trunk
http://mallet.cs.umass.edu/api/cc/mallet/util/Maths.html

*/


public class Kinship
	{
	private float treshold;
	private Model model;
	private Transaction transaction;
	private LinkageRow linkageRow;
	private Set<Integer> pair_ids=new TreeSet<Integer>();
	private SortedMap<Integer,IBDStatus> pairid2theoritical=new TreeMap<Integer, IBDStatus>();
	private SortedMap<Integer,IBDStatus> pairid2observed=new TreeMap<Integer, IBDStatus>();
	private SortedMap<Integer,PairOfIndividual> pairid2couple=new TreeMap<Integer, PairOfIndividual>();
	//private String rsName;
	private TreeMap<Float,Integer> typesAA = new TreeMap<Float, Integer>();
	private TreeMap<Float,Integer> typesAU = new TreeMap<Float, Integer>();
	private int nbPairsAA=0;
	private int nbPairsAU=0;
	public Kinship(Model model,Transaction txn,float treshold)
		{
		this.model=model;
		this.transaction=txn;
		this.treshold=treshold;
		}
	/**
	 * Copy model, transaction and treshold
	 * @param cp
	 */
	public Kinship(Kinship cp)
		{
		this.model=cp.model;
		this.transaction=cp.transaction;
		this.treshold=cp.treshold;
		}
	
	
	public float getTreshold()
		{
		return treshold;
		}
	
	
	public Model getModel()
		{
		return model;
		}
	
	public Transaction getTransaction()
		{
		return transaction;
		}
	
	
	
	private int getIbd0(int pair_id)
		{
		return this.linkageRow.get(pair_id).getIbd0() < getTreshold() ? 0:1;
		}
	
	
	public void setPairIds(Set<Integer> pair_ids)
		{
		this.pair_ids.clear();
		this.pairid2theoritical.clear();
		this.typesAA.clear();
		this.typesAU.clear();
		this.nbPairsAA=0;
		this.nbPairsAU=0;
		this.linkageRow=null;
		for(Integer pair_id:pair_ids)
			{
			IBDStatus status= getModel().getTheoriticalIBDStatusByPairId(getTransaction(), pair_id);
			if(status==null) throw new RuntimeException("No theo status for pair_id:"+pair_id);
			//discard those having a fixed probability
			if(status.getIbd0()==1f || status.getIbd1()==1f || status.getIbd2()==1f)
				{
				continue;
				}
			
			this.pairid2theoritical.put(pair_id,status);
			
			PairOfIndividual p= getModel().getPairOfIndividualByPairId(getTransaction(), pair_id);
			if(p==null) throw new RuntimeException("No pair of indis for pair_id:"+pair_id);
			this.pairid2couple.put(pair_id, p);
			this.pair_ids.add(pair_id);
			
			switch(p.getConcordance())
				{
				case AFFECTED:
					{
					Integer f= this.typesAA.get(status.getIbd0());
					this.typesAA.put(status.getIbd0(),1+(f==null?0:f));
					this.nbPairsAA++;
					break;
					}
				case DISCORDANT:
					{
					Integer f= this.typesAU.get(status.getIbd0());
					this.typesAU.put(status.getIbd0(),1+(f==null?0:f));
					this.nbPairsAU++;
					break;
					}
				}
			}
		}
	
	
	public Stat getByMarkerName(String rsName)
		{
		//this.rsName=rsName;
		this.linkageRow=getModel().getIBDByMarker(getTransaction(), rsName);
		this.pairid2observed.clear();
		for(Integer pair_id:this.pair_ids)
			{
			IBDStatus obs= this.linkageRow.get(pair_id);
			this.pairid2observed.put(pair_id, obs);
			}
		
		double zScoreAA=0;
		for(Float catIbd0: this.typesAA.keySet())
			{
			int eventsAA=0;
			int totalAA=0;
			for(Integer pair_id: this.pair_ids)
				{
				if(this.pairid2theoritical.get(pair_id).getIbd0()!=catIbd0) continue;
				
				if(this.pairid2couple.get(pair_id).getConcordance()!=PairOfIndividual.Concordance.AFFECTED) continue;
				totalAA++;
				
				if(getIbd0(pair_id)!=1) continue;
				eventsAA++;
				}
			
			double prAA = pbinom( eventsAA-1, totalAA, catIbd0,true,false ) +
					dbinom(eventsAA,totalAA, catIbd0,false)
					; 
			
			prAA=(prAA>=1.0?0.0999999:prAA);
			
			zScoreAA +=
					qnorm(prAA,0,1.0,false,false)*Math.sqrt(totalAA);
			}
		double saveZscoreAA=zScoreAA;
		zScoreAA=zScoreAA/Math.sqrt(nbPairsAA);//truc a utiliser LA ICI 
		
		
		
		
		double zScoreAU=0;
		for(Float catIbd0: this.typesAU.keySet())
			{
			int eventsAU=0;
			int totalAU=0;
			for(Integer pair_id: this.pair_ids)
				{
				if(this.pairid2theoritical.get(pair_id).getIbd0()!=catIbd0) continue;
				
				if(this.pairid2couple.get(pair_id).getConcordance()!=PairOfIndividual.Concordance.DISCORDANT) continue;
				totalAU++;
				
				if(getIbd0(pair_id)!=1) continue;
				eventsAU++;
				}
			
			double prAU = pbinom( eventsAU-1, totalAU, catIbd0,false,false) +
					dbinom(eventsAU,totalAU, catIbd0,false)
					; 
			
			prAU=(prAU>=1.0?0.0999999:prAU);
			
			zScoreAU +=
					qnorm(prAU,0,1.0,false,false)*Math.sqrt(totalAU);
			}
		double saveZscoreAU=zScoreAU;
		zScoreAU=zScoreAU/Math.sqrt(nbPairsAU);

		//
		
		
		Stat stat=new Stat();
		stat.zScoreAA=zScoreAA;
		stat.zScore= (saveZscoreAA+saveZscoreAU)/Math.sqrt(nbPairsAA+nbPairsAU);
		
		return stat;
		}
	
	
	  public static class Stat
	  	{
		private double zScoreAA=Double.NaN;
		private double zScore=Double.NaN;
		
		public boolean isNan()
			{
			return Double.isNaN(zScoreAA) || Double.isNaN(zScore);
			}
		
		public double lodScoreAA()
			{
			if(isNan()) return Double.NaN;
			return (this.zScoreAA*this.zScoreAA)/2.0*Math.log(10);
			}
		
		public double lodScore()
			{
			if(isNan()) return Double.NaN;
			return (this.zScore*this.zScore)/2.0*Math.log(10);
			}
	  	}
	
	  /**
	   *
	   * Calculates the value of the cumulative distribution function
	   *
	   * @param dist the distribution
	   * @param q the value
	   * @param lowerTail if true, return the value P(x < q), otherwise P(x > q)
	   * @param logP  if true, return the natural logarithm of the probability
	   * @return  the probability that the random variable will take the value less than (greater than)
	   * {@code q}
	   */
	  private static double _p(Distribution dist, double q, boolean lowerTail, boolean logP) {
	    double p;
	    try {
	      p = dist.cumulativeProbability(q);
	    } catch (MathException e) {
	      return Double.NaN;
	    } catch (MathRuntimeException e) {
	      return Double.NaN;
	    }
	    if (!lowerTail) {
	      p = 1.0 - p;
	    }
	    if (logP) {
	      p = Math.log(p);
	    }

	    return p;
	  }

	  private static double _d(IntegerDistribution dist, double x, boolean log) {
		    double d = dist.probability(x);
		    if (log) {
		      d = Math.log(d);
		    }
		    return d;
		  }


	  private static double _q(ContinuousDistribution dist, double p, boolean lowerTail, boolean logP) {
		    if (logP) {
		      p = Math.exp(p);
		    }
		    double q = 0;
		    try {
		      q = dist.inverseCumulativeProbability(p);
		    } catch (IllegalArgumentException e) {
		      return Double.NaN;
		    } catch (MathException e) {
		      return Double.NaN;
		    } catch (MathRuntimeException e) {
		      return Double.NaN;
		    }
		    if (!lowerTail) {
		      q = -q;
		    }
		    return q;
		  }


	
	private static double pbinom (int x, int size, double prob,boolean lowerTail,boolean logP) 
		{
		 return _p(new BinomialDistributionImpl(size, prob), x, lowerTail, logP);
		}
	
	private static double  dbinom(double x,  int size, double prob, boolean log)
		{
		return _d(new BinomialDistributionImpl(size, prob), x, log);
		}
		
	private static double qnorm(double p,double mean, double sd, boolean lowerTail, boolean logP)
		{
		return _q(new NormalDistributionImpl(mean, sd), p, lowerTail, logP);
		}
	}
